**********************
*****mata program*****
**********************  
mata
version 11
mata clear
mata set matastrict on

void notch_ct_roundq (
            string scalar    incbin_name,
            string scalar    freq_name, 
			real scalar      notch,       /* standardlized  */
            real scalar      low,     /* negative sign _number of bins below notch where excluded range for counterfactual starts */
            real scalar      high, 
            real scalar      degree,  /* degree of polynomial */
            real scalar      nboot,   /* number of iterations in bootstrap */
            real scalar      sf_s) 	

{
real matrix         Poly,   get_beta_dummy, pick_right,  pick_bunch_and_right , dummies, x_nodum,  x, 
                    xxinv,  xxinvx,   x_plot, rounds , get_beta_polyround,get_beta_dummy_bunch,get_beta_dummy_hole,sround,
					get_beta_polyroundsround , get_beta_poly,get_beta_round,  get_beta_sround ,x_poly,x_poly_round,
                    get_beta_roundsround
                       
        
real colvector    incbin, freq, round80, round20, round50, iota,  cons, poly1, poly2, agg,  beta_hat,  beta_hat_dum, y_counter,
                   y_new,  y_hat,  y_hat_pred, error_adj, b_boot, m_boot, bn_boot,  bf_boot,  draw_j, resample,  freq_j,
                   y_counter_nosround, y_counter_noroundsround, y_actual_nosround, y_actual_noroundsround ,
				   beta_hat_j,  beta_hat_dum_j ,  beta_hat_polyround, beta_hat_dummy_bunch,beta_hat_dummy_hole, 
                   beta_hat_dummy_hole_j,beta_hat_dummy_bunch_j,Dummy,Z100, W50,y_actual_noround,
				   beta_hat_poly, beta_hat_round, beta_hat_sround,beta_hat_polyroundsround,y_counter_noround ,
				   y_counter_j, y_counter_nosround_j,  y_counter_noroundsround_j,
				   beta_hat_polyroundsround_j,beta_hat_polyround_j,   beta_hat_poly_j,
				   y_roundsround_j ,  beta_hat_roundsround_j , x_roundsround, y_roundsround , beta_hat_roundsround
				   
				          
real scalar        mn, rows, se,  last_poly_b, last_polyround_b, first_dummy_b, ndum,  xcols, minim,  maxim, maxim_rx, notch_rx, 
                   first_bunch_rx,  last_bunch_rx,  last_left_rx, first_right_rx,  people_nearbunch, bn, bf, b, 
			       n_it, bn_j, bf_j, b_j, mp_j, m_j, n_it_j, i,j,ii, jj, b_se, m_se, nround, ndumbunch,ndumhole,first_hole_rx, 
				   last_hole_rx,mp,af, m, af_boot,af_j,nsround,
			       last_dzd_rx, num_dzd
				   
string scalar      b_s		
			  
/* GETTING MATRICES FOR ESTIMATION */	
  incbin			=st_data(.,incbin_name) 
  freq			    =st_data(.,freq_name)

  round50           =st_data(.,"round50")
  Z100              =st_data(.,"Z100")
  rounds            =round50, Z100
  
  
/* CONSTRUCTING CONSTANT, POLYNOMIAL TERM, AND DUMMIES TERM */
  rows            =rows(incbin)
  iota            =J(rows,1,1)
  cons            =iota			
			
/* NORMALIZE PLOY1 - DEMEAN AND DIVIDE BY SE   */
  mn              =(iota'*incbin)/(iota'*iota)     
  poly1           =incbin:-mn                  
  se              =(sqrt(cross(poly1,poly1)/rows)) 
  poly1           =poly1/se                     //normalized
  Poly            =poly1
  for (i=2; i<=degree; i++) Poly = Poly,poly1:^i	
  poly2           =poly1:^2
  sround=Z100:*poly1, Z100:*poly2, round50:*poly1,  round50:*poly2

/* ROW NUMBERS IN BETA */
  nround          =2
  nsround         =4
  
  ndumhole        =-low	        //# of hole dummies
  ndumbunch       =high+1
  last_polyround_b=1+degree+nround+nsround
  first_dummy_b   =1+degree+nround+nsround+1  
  
  /*defining notch and bunching window */	
  ndum  = high-low+1            //# of dummies
  num_dzd=high*2+1
  xcols =1+degree+nround+nsround+ndum 

/* ROW NUMBERS FOR BUNCHING RANGE */
  minim = colmin(incbin) 		 /* lowest incbin */
  maxim = colmax(incbin)        /* highest incbin */      


/* ROW NUMNBERS FOR EXCLUDED RANGE */	
  maxim_rx        = (maxim-minim)/sf_s+1
  notch_rx        =  notch-minim /sf_s+1     /*bin number of notch*/      
		
  first_hole_rx  = notch_rx+low	     /*bin number of first_bunch, usually 1 if minim=notch+low*/
  last_hole_rx   = notch_rx-1
  first_bunch_rx = notch_rx
  last_bunch_rx  = notch_rx+high		    
  last_left_rx   = first_hole_rx-1
  first_right_rx = last_bunch_rx+1
  last_dzd_rx    = notch_rx-high
  
  /* USEFUL SCALARS, VECTORS AND MATRICES */
   get_beta_polyroundsround      = diag(J((1+degree+nround+nsround ),1,1))  , J((1+degree+nround+nsround ),ndum,0)
   get_beta_dummy_hole           = J(ndumhole,(1+degree+nround+nsround),0) , diag(J(ndumhole,1,1))  , J(ndumhole,ndumbunch,0)
   get_beta_dummy_bunch          = J(ndumbunch,(1+degree+nround+nsround),0) ,J(ndumbunch,ndumhole,0), diag(J(ndumbunch,1,1))					 
   get_beta_roundsround          = J((nround+nsround),(1+degree),0)         , diag(J((nround+nsround),1,1)),  J((nround+nsround),ndum,0)  

  
  /* CREATING DUMMIES FOR EXCLUDED RANGE IN COUNTERFACTUAL */
  dummies = J(rows, ndum, .)
  ii=0
  for(jj=first_hole_rx; jj<=last_bunch_rx; jj++) {
       ii = ii+1
       dummies[.,ii] = e(jj,rows)'
    }
		
		
/* CREATING THE MATRICES NEEDED FOR ESTIMATION OF BETA_HAT */
  x_nodum = cons, Poly, rounds, sround
  x = x_nodum, dummies 
  xxinv = invsym(cross(x,x))
  xxinvx = xxinv*x'
  
  x_poly=cons, Poly
  x_poly_round=cons, Poly, rounds
  x_roundsround=rounds, sround

 /* CONSTRUCTING BETAHAT AND Y_COUNTER*/
  beta_hat                 =xxinvx*freq
  beta_hat_polyroundsround =get_beta_polyroundsround *beta_hat	
  beta_hat_dummy_hole      =get_beta_dummy_hole      *beta_hat
  beta_hat_dummy_bunch     =get_beta_dummy_bunch     *beta_hat
  beta_hat_roundsround     =get_beta_roundsround     *beta_hat

  
 /* CONSTRUCTING FITTED DATA */ 
  //ct     (no dum)
  y_counter               =x_nodum*beta_hat_polyroundsround		
  st_matrix("y_counter",y_counter)
  
  y_roundsround           =x_roundsround*beta_hat_roundsround 
   
  y_counter_noroundsround =y_counter -  y_roundsround
  st_matrix("y_counter_noroundsround",y_counter_noroundsround)
  
  //numbin2_rv2 (no sround, no round)
  y_actual_noroundsround  =freq      -  y_roundsround
  st_matrix("y_actual_noroundsround",y_actual_noroundsround)


 /* CALCULATING EXCESS MASS AND ELASTICITY/EVASION RESPONSE */
  mp = colsum(beta_hat_dummy_hole)
  bn = colsum(beta_hat_dummy_bunch)
  af=colsum(y_counter[(first_hole_rx)..(last_hole_rx)])/(-low)
  b = bn/(af)
  m = mp/af
  st_global("b", strofreal(b,"%9.2f"))
  st_global("m", strofreal(m,"%9.2f"))

 /* CONSTRUCTING NEW FITTED DATA */
  y_hat=x*beta_hat        //py
  
 /* BOOTSTRAPPING FOR STANDARD ERRORS */	  
  error_adj      = freq - y_hat    //res
  b_boot         = J(nboot,1,0)
  m_boot         = J(nboot,1,0)
	
 /* RUNNING THE BOOTSTRAP LOOP - GENERATE NEW DATASET */ 
  for (j=1; j<=nboot; j++) 		{
  draw_j = ceil(uniform(rows,1)*rows)    
  resample = error_adj[draw_j]
  freq_j = y_hat + resample
			
 /*BOOTSTRAPPING OTHER VARIABLES */ 
  beta_hat_j                   = xxinvx*freq_j
  beta_hat_polyroundsround_j   =get_beta_polyroundsround*beta_hat_j
  beta_hat_dummy_hole_j        =get_beta_dummy_hole*beta_hat_j
  beta_hat_dummy_bunch_j       =get_beta_dummy_bunch*beta_hat_j
  beta_hat_roundsround_j       =get_beta_roundsround*beta_hat_j

  
 /* CONSTRUCTING FITTED DATA */
  y_counter_j                  =x_nodum     *beta_hat_polyroundsround_j		
  y_roundsround_j              =x_roundsround*beta_hat_roundsround_j 
  y_counter_noroundsround_j    =y_counter_j -  y_roundsround_j

 
 /* CALCULATE EXCESS MASS MEASURES */
  bn_j = colsum(beta_hat_dummy_bunch_j)
  af_j = colsum(y_counter_j[(first_hole_rx)..(last_hole_rx)])/(-low)
  b_j  = bn_j/af_j	
  b_boot[j] =b_j
  
  mp_j = colsum(beta_hat_dummy_hole_j)
  m_j = mp_j/af_j
  m_boot[j] =m_j
   }  
   /* end bootstrap loop */

 /* WRITE RESULTS INTO STATA */
   b_se        =sqrt(quadvariance(b_boot))   
  st_global("b_se", strofreal(b_se,"%9.2f"))		

  m_se        =sqrt(quadvariance(m_boot))   
  st_global("m_se", strofreal(m_se,"%9.2f"))
  }
mata mosave notch_ct_roundq(),replace
end
